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Abstract 

A cycle expansion for the Lyapunov exponent of a product of ran- 
^ , dom matrices is derived. The formula is non-perturbative and numer- 

ically effective, which allows the Lyapunov exponent to be computed 
to high accuracy. In particular, the free energy and the heat capacity 
are computed for the one-dimensional Ising model with quenched dis- 
order. The formula is derived by using a Bernoulli dynamical system 
to mimic the randomness. 

5_j ■ The product of random matrices often appears in the study of disor- 
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dered materials and of dynamical systems. The physical quantities of these 
systems are related to the rate of growth of the random product — the Lya- 
punov exponent. For example, in the study of an Ising model with quenched 
randomness the Lyapunov exponent is proportional to the free energy per 
particle; in the Schrodinger equation with a random potential, the Lyapunov 
exponent is proportional to the localization length of the wave function; and 
in the motion of a classical particle, the Lyapunov exponent indicates the 
degree of sensitivity to initial conditions (chaos). Since Dyson |[| studied a 
system of harmonic oscillators with random couplings, many problems have 
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been reduced to the study of a Lyapunov exponent. In these problems the 
product of random matrices appears when a discrete version of a differential 
operator is considered, or when the problem is solved on a lattice. Further 
applications of Lyapunov exponents, and related derivations, are reviewed 
in the paper by Alexander et al. and the book by Crisanti, Paladin and 
Vulpiani Q. 

Few are the analytic results for the Lyapunov exponent of a product of 
random matrices, and very little has been determined about related systems 
without resorting to Monte Carlo simulations. A theorem of Oseledec [Q] 
states that the norm of the random product grows exponentially with the 
number of multiplied terms at a rate given by the Lyapunov exponent, but 
the theorem does not provide a method for determining the exponent. The 
two known methods for calculating the Lyapunov exponent, weak disorder 
expansions H ||, |?]] and microcanonical approximations ||, have limitations. 
The weak disorder expansion imposes conditions on the eigenvalues of the 
matrices and is difficult to carry out to high orders; and the microcanoni- 
cal method, while general, does not provide a systematic expansion and is 
difficult to apply to large matrices. In this Letter a formula for comput- 
ing the Lyapunov exponent will be derived. It is simple to evaluate and is 
non-perturbative in character, with the first few terms providing a good nu- 
merical approximation. In particular, it gives all thermodynamic quantities 
for the one-dimensional Ising model with a discrete valued random magnetic 
field when the disorder averaging is done over the free energy (this is the 
more difficult quenched disorder case). 

The formula is obtained by viewing the random product as a statistical 
mechanical system which is solved using the cycle expansion || of its ther- 



modynamical zeta function [10|. Cycle expansions have been very success- 



ful for obtaining non-perturbative expansions of chaotic dynamical systems 



[11, 12], of generalized Ising systems [13, 14 1, and of scattering problems in 
quantum mechanics Jlq ]. 

We will consider the product 

g w = n ^ 

0<fc<n 

of matrices chosen at random from a discrete set. This includes many 
cases of interest and can be used to approximate a continuous distribution. 



The maximal Lyapunov exponent 7 can be expressed [16] as the rate of 



exponential growth of the norm of the product with the number n of 
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matrices multiplied: 

7 = lim -{ln\\G {n) \\) (2) 

n — >oo fi 

The average is over all possible realizations of the product, each product 
taken with the appropriate probability. The theorem of Oseledec Q guar- 
antees that the limit exists for almost every realization. The definition of 
the Lyapunov exponent appears to depend on the matrix norm chosen, but 
it can be shown that its value remains unchanged as long as equivalent 
norms are used. For the finite dimensional vector space of n x n matrices all 
norms are equivalent, making the Lyapunov exponent independent of the 
norm chosen. 

Because it is very difficult to handle the logarithm inside the average 
we will substitute a power and a derivative for the logarithm, and write the 
Lyapunov exponent as 



(3) 

a=0 



7= lim d a -(\\G^\\ a ) 

n-^-oo 77, 

To determine the averages, introduce the generating function 

C(z,a)=exp(j2-(\\G {n) \\ a )) , (4) 
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which is the Ruelle zeta function 13] for a statistical mechanical system 
with (||G^ n ^|| a ) as the free energy. The zero z(a) of 1/C(z,a) gives the 
exponential growth of (||G (n) || Q ) (see Ref. pi Thm 5.29]), and by using the 
special values z(0) = 1 and c^C _1 (l,0) = —1, the Lyapunov exponent can 
be re-expressed as 

7 = -d a ln z(a = 0) = -daC 1 (1,0) . (5) 

The expression in terms of the derivative of the zeta function is of no 
advantage unless it can be computed in an efficient manner. If the terms 
of the zeta function satisfy certain combinatorial properties the inverse zeta 
function can be written as a cycle expansion ||, which is rapidly convergent 
and offers a practical scheme for evaluating the Lyapunov exponent. The 
average z n (\\G^ ||) is the sum of terms of the form 

t G = z n(G) Prob(G)||G|| Q (6) 

with n(G) being the number of matrices in the product G. If the weights 
are cyclic, as in tAAB = t-ABA = t-BAA, and multiplicative, as in tABAB = 
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(tAB) 2 > then the inverse zeta function can be expanded into a cycle expan- 
sion 

CHz,a)= Il(l-tc) (7) 

GeP 

with the product being over the set P of all possible prime products. A 
product of matrices is prime if it is not the repeat of a smaller length product. 
Two products are equivalent for the expansion if they differ by a cyclic 
rotation. For example, if AB is in the set P, then BA does not need to 
be in the set as it is equivalent to AB by cyclic rotation. Also ABAB and 
BAB A need not be in P as they are repeats of AB or of its cyclic rotation. 

To continue the derivation we use the independence of the Lyapunov 
exponent on the norm, and choose the most convenient norm for the cycle 
expansion. If we choose the norm absolute value of the largest eigenvalue 
(or eigenvalues), and write it in the peculiar form 

IIGir = lim \trG n \ a/n (8) 

n—*oo 

then it is simple to verify that the weights to are multiplicative and cyclic. 
They are multiplicative because 

/ \ <xp 

\\G p \\ a = lim \trG pn \ a/n = lim |trG pn | 1/(pn) = \\G\\ ap . (9) 

n^oo V pn^co J 

And they are cyclic because the trace is cyclic. The Prob(-) part of the 
weight is multiplicative and cyclic, as the product of numbers is. The cycle 
expansion for (||G^|| a ) is then obtained by expanding the infinite product 

C\z,a)= Yl (l-z n ( G )Prob(G)||G|| Q ) (10) 
GeP 

into a power series in z. This is essential when computing the zeros of C 1 - 
The power series in z converges as long as the matrices are hyperbolic (not 
all eigenvalues equal to one in modulus). The cycle expansion (|TJJ) can be 
used to compute the Lyapunov exponent in equation (||). In the case that 
there are two matrices forming the random product, A with probability p 
and B with probability q = 1 — p, the first few terms of the expansion of (||) 
are: 

-7 = pin + gin \\B\\ + pq (In \\AB\\ - In \\A\\ - In + 

+ ^(ln \\AAB\\ - In ||AB|| - In ||A||) + (11) 
+ pqq(\n \\ABB\\ - In \\AB\\ - In + ■ ■ • 
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method 



Lyapunov 



heat capacity 



Weak Disorder 
Micro canonical 
Monte Carlo 
Zeta (n = 15) 



1.5 

1.20 

1.1773 

1.17727361334268 



0.4 

0.3664 



Table 1: The Lyapunov exponent and its second derivative (proportional to 
the heat capacity) for the transfer matrices (see equation [12]) when j = 0.3 
and h = 1.4. The Monte Carlo and zeta function values are quoted so that 
errors are limited to the last digit. 



General expressions and a more detailed derivation can be found in Ref. 







One test of the expansion for the Lyapunov exponent is a product of ran- 
dom matrices that appears in the study of the one-dimensional Ising model 
with coupling constant J and a random magnetic field assuming the values 
±/i. The two matrices, chosen with equal probability and after factoring out 
a common term, are of the form 



1 ab k 
a b k 



(12) 



where k is +1 or —1, a is exp(— 2 J), and b is exp(— 2h). In this case both ma- 
trices have eigenvalues that are real and different from one, and are therefore 
hyperbolic. The Lyapunov exponent for these matrices can also be computed 
by Monte Carlo simulations, weak disorder expansions, and micro canonical 
approximations. Table |l| has the results for these methods. The Monte Carlo 
calculations corresponds to 128 realizations of a product of 10 6 matrices; the 
weak disorder expansion is carried to fifth order; there is no error estimate 
for the microcanonical method, except for a rigorous upper bound; and the 
value for the zeta function method was computed by including cycles up to 
length fifteen. The zeta function expansion takes 17 seconds on a Sparcsta- 
tion 1 computer, whereas if the Monte Carlo simulation had to be carried 
out to the precision obtained in the cycle expansion, it would require several 
hundred years of Sparcstation 1 time. The weak disorder expansion could in 
principle match the accuracy of the zeta function expansion if the terms of 
higher order were known. Also in table |] the second derivative of the Lya- 
punov exponent (proportional to the heat capacity at f3 = 1) is computed by 
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Figure 1: Number of digits that remain constant as longer cycles are included 
in the expansion of the Lyapunov exponent. The bullets (•) are for the 
random Ising model and the crosses (+) are for the degenerate 3x3 matrices. 
Also indicated in the plot are the accuracy of the weak disorder expansion 
(short dashes) and the microcanonical approximation (long dashes). 

numerical differentiation. Notice that, except for the zeta function, all ana- 
lytic methods fail to provide a value for the second derivative and that even 
though Monte Carlo does estimate the derivative with one digit, a better 
estimate would require a prohibitive amount of computer time. 

To study the convergence of the zeta function method one can plot the 
rate at which digits are gained in the value of the Lyapunov exponent as 
longer cycles are included in the expansion. The better a system is under- 
stood, the better the nature of the convergence. In figure [l] the number of 
correct digits as a function of the largest cycle considered is plotted. If 7„ 
is the approximation to the Lyapunov exponent when cycles up to length n 
are included, then the number of digits is defined as d(n) = log 10 (7 n _i — j n ). 
The straight line indicates that the convergence is exponentially fast in the 
length of the product. 

To further illustrate the method, figure || has a plot of the Lyapunov 
exponent for the Ising model pair of matrices chosen with equal probability 
and in units where J = h = (3. In the plot all points can be computed to 
machine precision, and the convergence rate is similar to that of figure |j. 
Thus the method is not limited to small values of the inverse temperature /3 




Figure 2: Lyapunov exponent (free energy) for the pair of matrices (12) 
that describe the Ising model with quenched randomness as a function of the 
inverse temperature (3. The dotted line is an interpolation of the computed 
points. 



as the weak disorder expansion, and can be used to obtain thermodynamic 
quantities at any temperature. 

The weak disorder expansion cannot be applied when there are repeated 
eigenvalues, so to illustrate the zeta function method the Lyapunov expo- 
nent has been computed for a pair of matrices with degenerate eigenvalues. 
The random products are formed from a pair of three by three matrices, 
which have the same eigenvalues, do not commute, and are not related by a 
similarity transformation. The eigenvalues are 2, 2, and 1. The exponential 
convergence of the method is not affected by the largest eigenvalue being de- 
generate, nor by the presence of an indifferent eigenvalue (the one of value 
one). In figure |l] the crosses are the plot of the number of non-changing 
digits of the Lyapunov exponent as a function of the cycle length. 

The cycle expansion developed in this Letter for the Lyapunov exponent 
is an efficient computational tool. It can be applied to a wide variety of 
matrix products without the limitations of other methods, and excludes only 
the matrices where all the eigenvalues are one in modulus. The method has 
been successfully applied to Ising models with a random magnetic field on 



a strip 17], and also in reproducing the branch point at zero temperature 
predicted by Derrida and Hilhorst [18] in the same model. 
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